{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import tensorflow as tf\n",
    "import datetime, os\n",
    "#hide tf logs \n",
    "os.environ['TF_CPP_MIN_LOG_LEVEL'] = '2'  # or any {'0', '1', '2'} \n",
    "#0 (default) shows all, 1 to filter out INFO logs, 2 to additionally filter out WARNING logs, and 3 to additionally filter out ERROR logs\n",
    "import scipy.optimize\n",
    "import scipy.io\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib.gridspec as gridspec\n",
    "from mpl_toolkits.axes_grid1 import make_axes_locatable\n",
    "from mpl_toolkits.mplot3d import Axes3D\n",
    "import time\n",
    "from pyDOE import lhs         #Latin Hypercube Sampling\n",
    "import seaborn as sns \n",
    "import codecs, json\n",
    "import matplotlib.transforms as mtrans\n",
    "\n",
    "# generates same random numbers each time\n",
    "np.random.seed(1234)\n",
    "tf.random.set_seed(1234)\n",
    "\n",
    "print(\"TensorFlow version: {}\".format(tf.__version__))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "88G3Lt8xn-Oo"
   },
   "source": [
    "# *Data Prep*\n",
    "\n",
    "Training and Testing data is prepared from the solution file"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "data = scipy.io.loadmat('Data/burgers_shock_mu_01_pi.mat')  \t# Load data from file\n",
    "x = data['x']                                   # 256 points between -1 and 1 [256x1]\n",
    "t = data['t']                                   # 100 time points between 0 and 1 [100x1] \n",
    "usol = data['usol']                             # solution of 256x100 grid points\n",
    "\n",
    "X, T = np.meshgrid(x,t)                         # makes 2 arrays X and T such that u(X[i],T[j])=usol[i][j] are a tuple"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "ZyGxyaOAcqpi"
   },
   "source": [
    "# *Test Data*\n",
    "\n",
    "We prepare the test data to compare against the solution produced by the PINN."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "yddknKA2Xohp"
   },
   "outputs": [],
   "source": [
    "''' X_u_test = [X[i],T[i]] [25600,2] for interpolation'''\n",
    "X_u_test = np.hstack((X.flatten()[:,None], T.flatten()[:,None]))\n",
    "\n",
    "# Domain bounds\n",
    "lb = X_u_test[0]  # [-1. 0.]\n",
    "ub = X_u_test[-1] # [1.  0.99]\n",
    "\n",
    "'''\n",
    "   Fortran Style ('F') flatten,stacked column wise!\n",
    "   u = [c1 \n",
    "        c2\n",
    "        .\n",
    "        .\n",
    "        cn]\n",
    "\n",
    "   u =  [25600x1] \n",
    "'''\n",
    "u = usol.flatten('F')[:,None] "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "aJ5oBRtEXnyu"
   },
   "source": [
    "# *Training Data*\n",
    "\n",
    "The boundary conditions serve as the test data for the PINN and the collocation points are generated using **Latin Hypercube Sampling**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "8UVJmvZbXjXb"
   },
   "outputs": [],
   "source": [
    "def trainingdata(N_u,N_f):\n",
    "\n",
    "    '''Boundary Conditions'''\n",
    "\n",
    "    #Initial Condition -1 =< x =<1 and t = 0  \n",
    "    leftedge_x = np.hstack((X[0,:][:,None], T[0,:][:,None])) #L1\n",
    "    leftedge_u = usol[:,0][:,None]\n",
    "\n",
    "    #Boundary Condition x = -1 and 0 =< t =<1\n",
    "    bottomedge_x = np.hstack((X[:,0][:,None], T[:,0][:,None])) #L2\n",
    "    bottomedge_u = usol[-1,:][:,None]\n",
    "\n",
    "    #Boundary Condition x = 1 and 0 =< t =<1\n",
    "    topedge_x = np.hstack((X[:,-1][:,None], T[:,0][:,None])) #L3\n",
    "    topedge_u = usol[0,:][:,None]\n",
    "\n",
    "    all_X_u_train = np.vstack([leftedge_x, bottomedge_x, topedge_x]) # X_u_train [456,2] (456 = 256(L1)+100(L2)+100(L3))\n",
    "    all_u_train = np.vstack([leftedge_u, bottomedge_u, topedge_u])   #corresponding u [456x1]\n",
    "\n",
    "    #choose random N_u points for training\n",
    "    idx = np.random.choice(all_X_u_train.shape[0], N_u, replace=False) \n",
    "\n",
    "    X_u_train = all_X_u_train[idx, :] #choose indices from  set 'idx' (x,t)\n",
    "    u_train = all_u_train[idx,:]      #choose corresponding u\n",
    "\n",
    "    '''Collocation Points'''\n",
    "\n",
    "    # Latin Hypercube sampling for collocation points \n",
    "    # N_f sets of tuples(x,t)\n",
    "    X_f_train = lb + (ub-lb)*lhs(2,N_f) \n",
    "    X_f_train = np.vstack((X_f_train, X_u_train)) # append training points to collocation points \n",
    "\n",
    "    return X_f_train, X_u_train, u_train \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "dp4nc2S7bwzz"
   },
   "source": [
    "# **PINN**\n",
    "\n",
    "Generate a **PINN** of L hidden layers, each with n neurons. \n",
    "\n",
    "Initialization: ***Xavier***\n",
    "\n",
    "Activation: *tanh (x)*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "class Sequentialmodel(tf.Module): \n",
    "    def __init__(self, layers, name=None):\n",
    "       \n",
    "        self.W = []  #Weights and biases\n",
    "        self.parameters = 0 #total number of parameters\n",
    "        \n",
    "        for i in range(len(layers)-1):\n",
    "            \n",
    "            input_dim = layers[i]\n",
    "            output_dim = layers[i+1]\n",
    "            \n",
    "            #Xavier standard deviation \n",
    "            std_dv = np.sqrt((2.0/(input_dim + output_dim)))\n",
    "\n",
    "            #weights = normal distribution * Xavier standard deviation + 0\n",
    "            w = tf.random.normal([input_dim, output_dim], dtype = 'float64') * std_dv\n",
    "                       \n",
    "            w = tf.Variable(w, trainable=True, name = 'w' + str(i+1))\n",
    "\n",
    "            b = tf.Variable(tf.cast(tf.zeros([output_dim]), dtype = 'float64'), trainable = True, name = 'b' + str(i+1))\n",
    "                    \n",
    "            self.W.append(w)\n",
    "            self.W.append(b)\n",
    "            \n",
    "            self.parameters +=  input_dim * output_dim + output_dim\n",
    "    \n",
    "    def evaluate(self,x):\n",
    "        \n",
    "        x = (x-lb)/(ub-lb)\n",
    "        \n",
    "        a = x\n",
    "        \n",
    "        for i in range(len(layers)-2):\n",
    "            \n",
    "            z = tf.add(tf.matmul(a, self.W[2*i]), self.W[2*i+1])\n",
    "            a = tf.nn.tanh(z)\n",
    "            \n",
    "        a = tf.add(tf.matmul(a, self.W[-2]), self.W[-1]) # For regression, no activation to last layer\n",
    "        \n",
    "        return a\n",
    "    \n",
    "    def get_weights(self):\n",
    "\n",
    "        parameters_1d = []  # [.... W_i,b_i.....  ] 1d array\n",
    "        \n",
    "        for i in range (len(layers)-1):\n",
    "            \n",
    "            w_1d = tf.reshape(self.W[2*i],[-1])   #flatten weights \n",
    "            b_1d = tf.reshape(self.W[2*i+1],[-1]) #flatten biases\n",
    "            \n",
    "            parameters_1d = tf.concat([parameters_1d, w_1d], 0) #concat weights \n",
    "            parameters_1d = tf.concat([parameters_1d, b_1d], 0) #concat biases\n",
    "        \n",
    "        return parameters_1d\n",
    "        \n",
    "    def set_weights(self,parameters):\n",
    "                \n",
    "        for i in range (len(layers)-1):\n",
    "\n",
    "            shape_w = tf.shape(self.W[2*i]).numpy() # shape of the weight tensor\n",
    "            size_w = tf.size(self.W[2*i]).numpy() #size of the weight tensor \n",
    "            \n",
    "            shape_b = tf.shape(self.W[2*i+1]).numpy() # shape of the bias tensor\n",
    "            size_b = tf.size(self.W[2*i+1]).numpy() #size of the bias tensor \n",
    "                        \n",
    "            pick_w = parameters[0:size_w] #pick the weights \n",
    "            self.W[2*i].assign(tf.reshape(pick_w,shape_w)) # assign  \n",
    "            parameters = np.delete(parameters,np.arange(size_w),0) #delete \n",
    "            \n",
    "            pick_b = parameters[0:size_b] #pick the biases \n",
    "            self.W[2*i+1].assign(tf.reshape(pick_b,shape_b)) # assign \n",
    "            parameters = np.delete(parameters,np.arange(size_b),0) #delete \n",
    "\n",
    "            \n",
    "    def loss_BC(self,x,y):\n",
    "\n",
    "        loss_u = tf.reduce_mean(tf.square(y-self.evaluate(x)))\n",
    "        return loss_u\n",
    "\n",
    "    def loss_PDE(self, x_to_train_f):\n",
    "    \n",
    "        g = tf.Variable(x_to_train_f, dtype = 'float64', trainable = False)\n",
    "    \n",
    "        nu = 0.01/np.pi\n",
    "\n",
    "        x_f = g[:,0:1]\n",
    "        t_f = g[:,1:2]\n",
    "\n",
    "        with tf.GradientTape(persistent=True) as tape:\n",
    "\n",
    "            tape.watch(x_f)\n",
    "            tape.watch(t_f)\n",
    "\n",
    "            g = tf.stack([x_f[:,0], t_f[:,0]], axis=1)   \n",
    "\n",
    "            z = self.evaluate(g)\n",
    "            u_x = tape.gradient(z,x_f)\n",
    "\n",
    "        u_t = tape.gradient(z,t_f)    \n",
    "        u_xx = tape.gradient(u_x, x_f)\n",
    "\n",
    "        del tape\n",
    "\n",
    "        f = u_t + (self.evaluate(g))*(u_x) - (nu)*u_xx\n",
    "\n",
    "        loss_f = tf.reduce_mean(tf.square(f))\n",
    "\n",
    "        return loss_f\n",
    "    \n",
    "    def loss(self,x,y,g):\n",
    "\n",
    "        loss_u = self.loss_BC(x,y)\n",
    "        loss_f = self.loss_PDE(g)\n",
    "\n",
    "        loss = loss_u + loss_f\n",
    "\n",
    "        return loss, loss_u, loss_f\n",
    "    \n",
    "    def optimizerfunc(self,parameters):\n",
    "        \n",
    "        self.set_weights(parameters)\n",
    "       \n",
    "        with tf.GradientTape() as tape:\n",
    "            tape.watch(self.trainable_variables)\n",
    "            \n",
    "            loss_val, loss_u, loss_f = self.loss(X_u_train, u_train, X_f_train)\n",
    "            \n",
    "        grads = tape.gradient(loss_val,self.trainable_variables)\n",
    "                \n",
    "        del tape\n",
    "        \n",
    "        grads_1d = [ ] #flatten grads \n",
    "        \n",
    "        for i in range (len(layers)-1):\n",
    "\n",
    "            grads_w_1d = tf.reshape(grads[2*i],[-1]) #flatten weights \n",
    "            grads_b_1d = tf.reshape(grads[2*i+1],[-1]) #flatten biases\n",
    "\n",
    "            grads_1d = tf.concat([grads_1d, grads_w_1d], 0) #concat grad_weights \n",
    "            grads_1d = tf.concat([grads_1d, grads_b_1d], 0) #concat grad_biases\n",
    "\n",
    "        return loss_val.numpy(), grads_1d.numpy()\n",
    "    \n",
    "    def optimizer_callback(self,parameters):\n",
    "               \n",
    "        loss_value, loss_u, loss_f = self.loss(X_u_train, u_train, X_f_train)\n",
    "        \n",
    "        u_pred = self.evaluate(X_u_test)\n",
    "        error_vec = np.linalg.norm((u-u_pred),2)/np.linalg.norm(u,2)\n",
    "        \n",
    "        tf.print(loss_value, loss_u, loss_f, error_vec)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "bOjuHdzAhib-"
   },
   "source": [
    "# *Solution Plot*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "UWqNuRMLhg4m"
   },
   "outputs": [],
   "source": [
    "def solutionplot(u_pred,X_u_train,u_train):\n",
    "    \n",
    "    fig, ax = plt.subplots()\n",
    "    ax.axis('off')\n",
    "\n",
    "    gs0 = gridspec.GridSpec(1, 2)\n",
    "    gs0.update(top=1-0.06, bottom=1-1/3, left=0.15, right=0.85, wspace=0)\n",
    "    ax = plt.subplot(gs0[:, :])\n",
    "\n",
    "    h = ax.imshow(u_pred, interpolation='nearest', cmap='rainbow', \n",
    "                extent=[T.min(), T.max(), X.min(), X.max()], \n",
    "                origin='lower', aspect='auto')\n",
    "    divider = make_axes_locatable(ax)\n",
    "    cax = divider.append_axes(\"right\", size=\"5%\", pad=0.05)\n",
    "    fig.colorbar(h, cax=cax)\n",
    "    \n",
    "    ax.plot(X_u_train[:,1], X_u_train[:,0], 'kx', label = 'Data (%d points)' % (u_train.shape[0]), markersize = 4, clip_on = False)\n",
    "\n",
    "    line = np.linspace(x.min(), x.max(), 2)[:,None]\n",
    "    ax.plot(t[25]*np.ones((2,1)), line, 'w-', linewidth = 1)\n",
    "    ax.plot(t[50]*np.ones((2,1)), line, 'w-', linewidth = 1)\n",
    "    ax.plot(t[75]*np.ones((2,1)), line, 'w-', linewidth = 1)    \n",
    "\n",
    "    ax.set_xlabel('$t$')\n",
    "    ax.set_ylabel('$x$')\n",
    "    ax.legend(frameon=False, loc = 'best')\n",
    "    ax.set_title('$u(x,t)$', fontsize = 10)\n",
    "    \n",
    "    ''' \n",
    "    Slices of the solution at points t = 0.25, t = 0.50 and t = 0.75\n",
    "    '''\n",
    "    \n",
    "    ####### Row 1: u(t,x) slices ##################\n",
    "    gs1 = gridspec.GridSpec(1, 3)\n",
    "    gs1.update(top=1-1/3, bottom=0, left=0.1, right=0.9, wspace=0.5)\n",
    "\n",
    "    ax = plt.subplot(gs1[0, 0])\n",
    "    ax.plot(x,usol.T[25,:], 'b-', linewidth = 2, label = 'Exact')       \n",
    "    ax.plot(x,u_pred.T[25,:], 'r--', linewidth = 2, label = 'Prediction')\n",
    "    ax.set_xlabel('$x$')\n",
    "    ax.set_ylabel('$u(x,t)$')    \n",
    "    ax.set_title('$t = 0.25s$', fontsize = 10)\n",
    "    ax.axis('square')\n",
    "    ax.set_xlim([-1.1,1.1])\n",
    "    ax.set_ylim([-1.1,1.1])\n",
    "\n",
    "    ax = plt.subplot(gs1[0, 1])\n",
    "    ax.plot(x,usol.T[50,:], 'b-', linewidth = 2, label = 'Exact')       \n",
    "    ax.plot(x,u_pred.T[50,:], 'r--', linewidth = 2, label = 'Prediction')\n",
    "    ax.set_xlabel('$x$')\n",
    "    ax.set_ylabel('$u(x,t)$')\n",
    "    ax.axis('square')\n",
    "    ax.set_xlim([-1.1,1.1])\n",
    "    ax.set_ylim([-1.1,1.1])\n",
    "    ax.set_title('$t = 0.50s$', fontsize = 10)\n",
    "    ax.legend(loc='upper center', bbox_to_anchor=(0.5, -0.35), ncol=5, frameon=False)\n",
    "\n",
    "    ax = plt.subplot(gs1[0, 2])\n",
    "    ax.plot(x,usol.T[75,:], 'b-', linewidth = 2, label = 'Exact')       \n",
    "    ax.plot(x,u_pred.T[75,:], 'r--', linewidth = 2, label = 'Prediction')\n",
    "    ax.set_xlabel('$x$')\n",
    "    ax.set_ylabel('$u(x,t)$')\n",
    "    ax.axis('square')\n",
    "    ax.set_xlim([-1.1,1.1])\n",
    "    ax.set_ylim([-1.1,1.1])    \n",
    "    ax.set_title('$t = 0.75s$', fontsize = 10)\n",
    "    \n",
    "    plt.savefig('Burgers.png',dpi = 500)   "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "YRuuEXx-eeWa"
   },
   "source": [
    "# *Model Training and Testing*\n",
    "\n",
    "A function '**model**' is defined to generate a NN as per the input set of hyperparameters, which is then trained and tested. The L2 Norm of the solution error is returned as a comparison metric"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "N_u = [20,40,80,100,200] #Total number of data points for 'u'\n",
    "N_f = [2000,4000,8000,10000] #Total number of collocation points \n",
    "\n",
    "runs = 5\n",
    "\n",
    "#store errors\n",
    "data = np.zeros((runs,len(N_f),len(N_u))) # depth(10) x rows (4) x columns(5)\n",
    "\n",
    "for k in range(runs):\n",
    "    for i in range(len(N_f)):\n",
    "        for j in range(len(N_u)):\n",
    "            \n",
    "            # Training data\n",
    "            X_f_train, X_u_train, u_train = trainingdata(N_u[j],N_f[i])\n",
    "\n",
    "            layers = np.array([2,20,20,20,20,20,20,20,20,1]) #8 hidden layers\n",
    "\n",
    "            PINN = Sequentialmodel(layers)\n",
    "\n",
    "            init_params = PINN.get_weights().numpy()\n",
    "\n",
    "            start_time = time.time() \n",
    "\n",
    "            # train the model with Scipy L-BFGS optimizer\n",
    "            results = scipy.optimize.minimize(fun = PINN.optimizerfunc, \n",
    "                                              x0 = init_params, \n",
    "                                              args=(), \n",
    "                                              method='L-BFGS-B', \n",
    "                                              jac= True,        # If jac is True, fun is assumed to return the gradient along with the objective function\n",
    "                                              callback = None, \n",
    "                                              options = {'disp': None,\n",
    "                                                        'maxcor': 200, \n",
    "                                                        'ftol': 1 * np.finfo(float).eps,  #The iteration stops when (f^k - f^{k+1})/max{|f^k|,|f^{k+1}|,1} <= ftol\n",
    "                                                        'gtol': 5e-10, \n",
    "                                                        'maxfun':  50000, \n",
    "                                                        'maxiter': 1,\n",
    "                                                        'iprint': -1,   #print update every 50 iterations\n",
    "                                                        'maxls': 50})\n",
    "\n",
    "            elapsed = time.time() - start_time                \n",
    "            #print('Training time: %.2f' % (elapsed))\n",
    "\n",
    "            PINN.set_weights(results.x)\n",
    "\n",
    "            ''' Model Accuracy ''' \n",
    "            u_pred = PINN.evaluate(X_u_test)\n",
    "\n",
    "            error_vec = np.linalg.norm((u-u_pred),2)/np.linalg.norm(u,2)        # Relative L2 Norm of the error (Vector)\n",
    "            data[k][i][j] = error_vec\n",
    "            \n",
    "            print('Test Error: %.5f'  % (error_vec))\n",
    "            \n",
    "        print('-----') #print at the end of inner-most loop \n",
    "                "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 84,
   "metadata": {},
   "outputs": [],
   "source": [
    "#convert 3d array to 2d\n",
    "data_2d = np.reshape(data, (runs*len(N_f),len(N_u)))\n",
    "\n",
    "np.savetxt('trend.txt', data_2d)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "data_2d = np.loadtxt('trend.txt')\n",
    "\n",
    "#convert 2d array to 3d\n",
    "data_3d = np.reshape(data_2d, (runs,len(N_f),len(N_u)))\n",
    "\n",
    "mean_2d = np.mean(data_3d, axis=0)\n",
    "\n",
    "std_2d = np.std(data_3d, axis=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "N_u = [20,40,80,100,200] #Total number of data points for 'u'\n",
    "N_f = [2000,4000,8000,10000] #Total number of collocation points \n",
    "\n",
    "runs = 5"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 88,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYwAAAEHCAYAAAC9TnFRAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAABF0UlEQVR4nO3dd3hUVfrA8e/JpCcQSiCUgAlVEaQ3aUE6S1FEBAvLiuJaFxYXRUVZsICyigqiSNtVEBWliAhICTXSg+IPEASEYOhSQgtJzu+PMwkpk2SSTMlM3s/z3IfMbfPOzXDfnHua0lojhBBC5MfH3QEIIYTwDJIwhBBC2EUShhBCCLtIwhBCCGEXSRhCCCHs4uvuAJxBKdUb6F2qVKnH6tSpU6hzXL58mZCQEMcGJnF4fAwSh8RR3GNwRBw7duw4o7WukGOD1tprl6ZNm+rCWrt2baGPdSSJo3jFoLXEkZ3EUbxi0LrocQDbtY17qjySEkIIYRdJGEIIIewiCUMIIYRdJGEIIYSwiyQMIYQQdpGEIYQQwi6SMIQQQthFEoYQQgi7SMIQQghhF0kYQggh7OIxY0kppUKAD4FkIFZrPdfNIQkhRIni1hKGUmqWUuqUUmpPtvXdlVL7lVIHlVIvWFf3AxZorR8D+rg8WCGEKOHc/UhqDtA98wqllAWYCvQA6gGDlFL1gEjgmHW3VBfGKIQQAjcnDK31euBcttUtgINa60Na62RgPtAXSMAkDXB/ohNCiBJHmZFs3RiAUlHAUq11fevr/kB3rfWj1tcPAy2B54EpwDVgY251GEqpYcAwgIiIiKbz588vVFxJSUmEhoYW6lhHkjiKVwwSh8RR3GNwRBwdO3bcobVulmODrTHPXbkAUcCeTK/7AzMyvX4YmFKYcxd6PowOHfSfDRsW7lgH85bx9b0lBq0ljuwkjuIVg9Ylaz6M40C1TK8jrevsppTqrZSafuHCBYcGJoQQJVlxTBjbgNpKqWillD8wEFhSkBNorb/VWg8LCwtzSoBCCFESubtZ7edAHFBXKZWglBqqtU4BngZWAHuBL7XWv7gzTiGEEG7uuKe1HpTL+mXAssKeVynVG+hdq1atwp5CCCFENsXxkVSRySMpIYRwPI8ZGsSl/ownNCXF3VEIIUSx4pUljKK2krp85RIcvgwnTjg4MiGE8FxemTCK+kjK/2QalqvA+PGODUwIITyYVyaMIjl6FN9LoABmzZJShhBCWHllwijSI6kJE27+fP06vPKK4wITQggP5pUJo9CPpBITYfZsU7owJ4KZMyEhwdEhCiGEx/HKhFFo48dDWlrWdWlp0KWLSR5CCFGCScLILC4OkpNzrt+3D15+2fXxCCFEMeKVCaPQdRi7dpF48Q/8XwE1FoJfC+LExT9g2DB44w344AOnxCuEEJ7AKxNGUZrV/nvdv7lhrcRITk3m3+vGwYcfwj33wD/+AV9+6eBohRDCM3hlwiisxEuJ/Hf3f0mv9U7VqUzfOZ340z/DvHnQti089BCsXu3eQIUQwg0kYWQyfv140nQalSwQWxUiLJCm02g9ozUrEtbBkiVQt64pbeza5e5whRDCpSRhZBKXEEdyajJjykHbIBhT7ua27nO7M3r7RG589y2ULQs9esChQ+4LVgghXMwrE0ZhK713Pb4L/dzvPBEGFgVPlQ9Cj0rk7PNneazJY0zYNIGYVQ/yx4LZkJICXbvCyZNO+hRCCFG8eGXCKNJYUr9MgPSuGKnX4edxBPsFM733dOb1m8dPJ3+i/pr+bJj6PPzxB/TsCZcuOTR+IYQojrwyYRTa1UQ4PBtlSV+RBr/NgKtmPKlBDQaxc9hOospE0X7fKD4Z3Q29ezf062e7/4YQQngRSRiZ/TwedLae3voG/PhIxsva5WsTNzSOZ1o8w7C0Rfz74WqwahUMGZKzl3gRlRkeQ68F0mFQCFE8SMLI7GwcpNkoKSQuhz93Z7wM8A3g/R7v8/WAr3nv1vO82j0APv8cRo6UIUSEEF5LZtzLrIdpKnuhji8aKPNrClz5A1a0gHW9oNtWCKqcsXu/2/rRpHITBpa7nzLntzJi8mRuRITj98JLbvoAQgjhPF5ZwijqjHtZBFeBDt/C9XOwri+kXMmyOapMFBse2ciJsc8xrz74jX6ZxCkTcjmZKKyYGBg+vJG7wxCiRPPKhFHUGfdyKNcY2syDc9sh7q856jn8LH5M7PY2ZecvIraWHxX+MZo1U59zzHsLIUQx4ZUJwyki+0Ljt+DYAvjpVZu79Li9L7Vjf+JgtVBajvgPr73Vi8vJlwv9lksXxPP9wp8LfbwQQjiS1GHY8PM/Q9Fa0y77hltHwsV98MtrULouRD+U49iqVW8l5cf9XGhanyfGfseAS3cw8ZnF1K9Y3yWxCyFKtpgYOH++EfHxjj+3lDBsuPrHXfzj7dic03krBc0+hIoxsGUonN5k83jfSlUov34boSFl+XjKEfr+pxkzds5ASwsqIYQHk4RhwyerBhN/pCHjx9vYaPGHdl9DyC2w/m5IymU8qZo1CVi5miqpwaya58tzXz7Gg988yMXrF50ZunAyqXwXJZkkjGwSE2Hh1r+gtQ8zZ8KBAzZ2CigHHZaCToV1vSE5l9ZYjRvjs2gxUSeT+en7W1gcP5+m05uyK1FGuhVCeB5JGNmMHw8pqaZq5/p1qFMHwsOhRQsYOBBefBFmzIA12+uQWOtr9MVfYeMASEuxfcK77kJ99hnVfz7K79vacv36FVrNbMXUrVPlEZUQwqNIpXcmiYkwezZkzKAE+PnBX/4CJ07Ajh3w9ddmoFqjI4/dNY3pQx9j8avDWXZqCjVqQI0aEB1t/i1bFtSAAXDyJOHPPsv+qoO5764zPP3906w5soaZfWZSJrCM6z+sEEIUkFcmDKVUb6B3rVq1CnTc+PE5h4NSCkJDYcUK8zo1FRISzFQYhw/DoUOPsvS3/fS9bRJxv9zKC9OfznJ8WBjWJPIMT7ZM5K5ZbzL5+ivUH3AX/9n9Ao0/bsz8e+fTMrJljnjOpEYw9txMlp+ASpUK9FGEEMLhvDJhaK2/Bb5t1qzZYwU5Li4u56CzycmwefPN1xYL3HKLWTp2tK5MmwAbfmXCff/glbdrcfBydw4dypxU4JdfoOeh15nKSYbOHcf5udNIrbaRYwMG0vqTtjS78Ca9w/9JrZo+GSWUORf/xU/JrRk/HqZOLdo18XTXr8PBg6GckOQphNt4ZcIorPRZV+MnxZCSkkKzFzbad6CPBe6cC6vaEbxjAHd03cwdd+Tsd5GWpkg89jFnHzjFh3FPEdPqKxad2cWqK0PZFvEvtu2NhTfnwJVw6xFDAcW0afDttxAYCL6+Nxc/v6yvC7KuoMf/+mtFzpxx7Hv5+poSnD12Xl5E8uU+kjyFcCNJGI7iF2rGnFrRwrSc6rYFAitm2cXHB6re4gs/fAGdOzNwyQMMXLEC3X4B07ZPY4RlBGUbNmJEtc/5dmo7Nu9OQPd/CL7+gqCgSjRtCjdumDqU7Mv163D5ctZ1tvbNvu7GDXs/YD2HXzIwJbb8EgtA8sE+gA8ffWQeG95+O0RF3VxCQ50SnhAex5mlcUkYjhQcCe2XwKr2po9GpzVgCbSxX7ApMrRrB337otav58nmT9I6sjUDFgxg9M6BWLYdQXd5E6pvRLcbz7G1U1m3zjmPY9LS8k8umzdvpUmTFrkmLHuSU27r8tt3y5assX7yialLyiw8/GbyiI7OmkyioswlF8KTaW2+/7b+/2RefvsNLl+2MG4cfPihY2OQhOFo5ZtB6//Bxvvgx6Fw52e2n7uUL29q0lu3hu7dYfNmGkc1ZuewnTTp8yMHgxKg0WzwSYPGs0mJG8P48ZWc8jjGxwcCAsySm4SEK9R3w+gmiYmm0UDmFuD+/vDjj3DlChw5YuqJjhwxy08/wZIlOeuiKlbMPaHccgsEBdkXj9SluE/mm2Vqas4bZWJiIAcO5H9DzW+xdW57l+PHb+Ojj5x3fvsp5syBV15x7PdUEoYzVO8PDV+H3S9B6VuhwRjb+1WrZpJG27bQrRts2kSp8HBCT3WGtk+BsjbZUqncaD2eb1eN4/kLl6keVt11n8XNbLVcS02Fjz82dRmtWuU8Ji3NNINOTyKZE8rOnbBwYc5HcZUq5Z1Q0pPp77+bv95cXZeidc4bzYULfpw44fqbYvbl2LE6/Pe/zjt/+pJ/tyUbXwYHy/4INfty40ZpSpWyvc1iMd+jvI4v7JIe19y5sGaNiTU1FYd/T5U3dx5r1qyZ3r59e4GPK3Clty1am6HQj3wKbebDLffnvu/GjdClC9xxB6xZQ2LaRWq8X4NrKdds7l6nfB261OhClxpd6BjdkdIBpQsfp51iY2OJiYlx+vtk17gxNgdRa9ToZiOFgkpNNSWX9CSSPakcPZrzr7kqVcyS/nWyWODhh01px9E3Rls33OyP4NzBx8f2zSot7TrBwQE5bl7OvDHaWg4c2EuDBrc55dzp2/NrpOGu/ydwszR+LdNtIyjItNIsaClDKbVDa90s+3opYTiLUtDyE7h8GH4cAiFREJ6zrwVgShhffAH33AP33cf4J6qTlm3ODX+LP/fceg8tq7bkh0M/MDt+NlO3TcWiLLSMbJmRQFpUbYGfxc/pH89V0pNCmeEmiSdNKUISt7JYIDLSLG3b5tyemgrHj+dMKOl9cdL3mT8fypSx/4YUHOzYG9jhwwe47bbaTrk52trfJ5dxIWJj49x2k8wax0liYm5zdxhuk1tp3JGlDEkYzmQJgHYLYWVLWN/XtJwKucX2vn36mOcsjz1GXKOyJAdkfQifnJrM/rP7md9/PiNaj+B6ynXiEuL44bcf+OHQD4xbN45/r/s3pfxL0TG6Y0YCqVO+DsretqsCMDfH6tXN0r69WZeYaBJEZkqZhOauuozY2OPExNR2z5uLYseefmRFJQnD2QLDzUCFK1ub5rZdNoFfKdv7PvoonDjBrjFjYNQoLn41ieDENHwPJ+a4KwX4BhATFUNMVAyvd3qdc1fPsebwmowEsmT/EgCqla5mkkfNLnSu0Znw4HBnf2Kv5Iq/3oQoivTSuJkP4zzx8WUc/h6SMFwh7DZo+yXE9oRNg6D9YtPZz5aXXjI1tm+9RXAAWK5j112pXFA5+tfrT/96/QH47dxv/HDIJI9v9n3DrPhZADSp3CSj9NGmehsCfW00+xU5uOKvNyGKO48ZrVYpVUMpNVMptcDdsRRK5a7Q7AP44zvY9a/c91MK3nsP/vIXfK9bh0GcPZucsznlrWa5mvy92d/5esDXnP7XaeKGxjEuZhwhfiH8J+4/dP60M2UnlqXbZ92YtHkSu0/sltFz87Brl2nH0KEDNGx4Hq3N68JWvAvhiVxSwlBKzQJ6Aae01vUzre8OvAdYgBla6wm5nUNrfQgY6rEJA6D2E3BhH+x/10zxWvtx2/tZLFClChprwkhOLtKzD18fX1pFtqJVZCvGdBjDpeuXWPf7uozHV//6wSSwiiEV6Vyjc0YJpGrpqoV6PyGEd3LVI6k5wBTgf+krlFIWYCrQBUgAtimllmCSx5vZjn9Ea33KNaE6WZN3IOkgbH8KStWESp1z7pOYCJ9+enOQ9dRUmD4dXn4ZKlcucgilAkrRq04vetXpBUDCxQRWHVrFD4d+YNWhVcz7eR4At4XfRteaXelSows+qR5TGBVCOIlLEobWer1SKirb6hbAQWvJAaXUfKCv1vpNTGmkUJRSw4BhABEREcTGxhb4HKEpKWitC3WsPSxpT9LYspfAtXezM/xDrvhl7YhX+913qZySkvV5YUoKl9q2Zccnn+TevrEIoojisXKPMbTsUA5dPsSOP3ew/c/tTNs2jfe2vIev8qXez/VoVrYZTcs2pW6pulhULvUwTpCSkoLGeb8Te50/34jU1FS3xwGQlJQkcRSzOIpDDE79jmqtXbIAUcCeTK/7Yx5Dpb9+GJiSx/HlgY+A34DR9rxn06ZNdWHseruD3vZmm0Ida7dLh7X+uqLWi2tqffV01m2NGqU/Is+5DBmi9Y0bzo0tk6s3rupVv63Sg2YO0k0+bqIZi2YsusyEMrrfF/30tG3T9MGzB50eR9g/OuiQp5z8O7FDhw5aN2z4p7vD0FprvXbtWneHoLWWOIpbDI74jgLbtY17qse0ktJanwX+7u44HCY0CtotgtUdYUM/uOsH028DMmpSL9TxRQNlfrWOizBuHIwdCxcvwrx5eQ/+5CCBvoF0qtEJy1ELMTExnLlyhtWHVt9sgbX3GwCiy0RnNN+9K/ouygWVc3psQgjXcmfCOA5Uy/Q60rquyAo7457LVWgNrWbD5gdg6+Pm59w62SkFr75quhYPHw69e5tBkUJCXBkx4cHh3F//fu6vfz9aa349+2tG8vh8z+dM3zkdH+VD08pNMxLIndXuxN/i79I4hRCO586EsQ2orZSKxiSKgcADjjixLuSMe24RNQgu/Qo/jzX9Neo9n/f+//iHmfd16FAz/tR335mJw91AKUXd8LrUDa/L0y2e5kbqDbYe35qRQCZumsgbG98gxC+EDlEdMlpf1atQT3qfC+GBXNWs9nMgBghXSiUAr2qtZyqlngZWYFpGzdJa/+Kg9/OMEka6+q/Axf0Q/wKUqg3V+uW9/5AhULo0DBpkunWuWFEsxtr2s/jRpnob2lRvw9iYsVy4doHYI7EZCWTZgWUAVClVJaP5bucanakU6v7YhRD5c1UrqUG5rF8GLHPC+3lOCQPM46ZWsyDpMGx+CLpsgHJN8z6mXz9Tuujb10zEtGqVGYe7GAkLDKPvrX3pe2tfAI5eOJrR9+O7X7/jf7tNK+sGFRtkPL5qf0t7gv1ktiMhiiNpXF9cWAKh/SIIqADr+sAVO6pzOnc2ieLMGTPs6r59Tg+zKKqHVWdok6HM7z+fU/86xY5hO5jQaQIVQyoyddtUesztQdmJZbnrv3cxYeMEdvyxI2PU3jR1nathP3MiqWA93oUQjuOVCUMp1VspNf3ChQvuDqVggiIgZincuATreqPKphE6JBWu5nGTbN0a1q0zMwK1a2dmCPIAPsqHJpWb8Hzb51k1eBXnnj/HiodW8GyLZzl39RyjV4+m2SfNqPh2Re5fcD9XQvaS5neR8evGuzt0IUosr0wYWutvtdbDwsLC3B1KwZVpYCZcOr+bkIc1lurAnnxuknfcARs2mBZTHTuanz1MsF8wXWt25e2ubxP/93hOjDzBZ/d8Rq86vVh3ZB2pfhdAwYxdM6SUIYSbeGXC8HhVe0KDsVjKgfIBDs3Ou5QBULu2mbmvShXo2hWWObxqyKUiQiN48I4HmXP3HO659R7QplVVcmoyz618zs3RCVEyeWXC8NhHUpld+QOdPi1n6tW8R7hNFxkJ69dDvXqmMvyLL5waoiskXkpkzu45oG6OpDvv53kcPHvQfUEJUUJ5ZcLw6EdSAFcT4fAcsgzVdOQzOPx5/sdWqGBmgb/zTtPsdvp0p4XpCuPXj88xXa1G0+WzLjnWCyGcyysThsf7eTzYuhnGPQB7Xre9LbOwMFi+HHr0gMcfh7feck6cLhCXEEdyanKO9UfOH2HcunFuiEiIkksSRnF0Ng7Sct4k8SsLP71sxp5KzudxW1CQGTpk4EB4/nl48UUzHpWH2fX4LvSrmrA/OxByug36VU3aK2kMaTSEf6/7d8ZYVkII5/OYwQcLwuN6emfXw8bgg2Bu+L9+ADtHwooW0H4hhNXL/Tz+/vDZZ6bE8eabcP48TJnilOHRXUkpxUd/+Yh9Z/YxeOFgaperTYOIBu4OSwiv59l3jlwUtQ5j+NJYHp2/1MFROYBSUPdZ6LQablwwSeNoPhMQWiwwbZopZUybBg8/bPpseLgA3wC+GfANYYFh9JnfhzNXzrg7JCG8nlcmDK9XsT103wFhDWDjfbBrFKSl5L6/UjBhgillzJtnhhW5etV18TpJ5VKVWXj/QhIvJTLgqwHcSPX8RChEcSYJw1MFV4XO66D2k7D3bVjbDa6dzvuYF16ADz80Y1D16GHm1fBwLaq2YHrv6aw9spaRK0e6Oxwh3C42FiZPjnfKuSVheDKLPzSfaubROL0JljeFs9vyPuaJJ0y9xsaN0KmTGYfKww1uOJh/tvonH2z9gJk7Z7o7HCG8llcmDK/ouFcQNYZA182mW/gPbeG3fG6aDzwAixbBnj3QoQMcd8i8VW41sctEutTowhPfPcHmY5vdHY4QXskrE4bHd9wrjHJNTL1GxQ6w5VEzg1/q9dz379XL9NU4dsyMdPvbb66L1Ql8fXyZ338+1cOq0++LfiRcTHB3SEJ4Ha9MGCVWQHmI+R7qjYaD02FVe7iSx42zQwfTK/zSJZM09uxxXawF1KgR1KqVlOc+5YLKsWTQEi7fuMw9X9zD1RueX7EvRHEiCcPb+Fig0RvQ7mu48H/wfRM4GZv7/s2amfGnfHygfXvYssVloTpDvQr1mNtvLtv/2M6wpcPQHthZUYjiShKGt6rWD7ptM6WONZ1h7zu59/SuV89UgpcrZyrCV692bawO1qduH8Z3HM9nP33GO3HvuDscIbyGJAxvFnYrdNsKkX1h10jYNAhSLtveNzrazKMRHQ09e5pKcQ/2UruXuK/efYxaNYoVB1e4OxwhvIIkDG/nVwraLoBGE+DYV7CiFVw8YHvfypXN7H2NG0P//vDpp66N1YGUUszuO5v6Fetz/4L7+fXsr+4OSQiP55UJo6jNap3Z8cUtlIJ6z0PMcriWCCuaw/Fchj4pV87MEx4TA4MHwwcfuDRURwrxD2HxwMX4+vjSd35fLl4vekdFr/tuCFEAXpkwSmSzWntU7mKa3obWhHW94adXbQ+VHhoKS5fC3XfDs8/Ca6955Ei3AFFlolgwYAEHzx3kwW8elDk0hCgCr0wYIg8ht0CXjaaz355xJnEk/5lzv8BA+OorU8oYM4aa06Z5bNKIiYrhve7vsfTXpYxZM8bd4QjhsbxyeHORD98gaDkLyreEHc/C8mbQbiGUvSPbfr4wezaEhVHtgw+gVCkzg5/FYvu8xdgTzZ4g/kQ8b2x8gzsi7uD++ve7OyQhPI6UMEoqpaD236Hzeki9BitbwZF5Offz8YH33uPI4MEwa5aZkOl6Hj3IiymlFFN6TqFNtTb8bfHfiD8R7+6QhPA4kjBKuvBWpl6jfHPY/CDsGA5p2YYJV4ojf/sbvPMOLFgAffrA5Vya5xZj/hZ/vh7wNeWDy9N3fl9OXT7l7pCE8CiSMAQEVYK7VkHd4bD/PVjdCa6eyLnfiBEwc6ZpRdW1q5nBz8NEhEaw6P5FnLp8iv5f9rc5X7gQwjZJGMLw8YOm78Kd8+DcdjNU+um4nPs98gh88QVs2wYdO8LJk66PtYiaVmnKzD4z2XB0A8OXD3d3OEJ4DLsThlLqS6XUV0qpt5VSg5RSdZ0ZmHCTqEHQ9UewBMHqDnDARuuo/v3h229h/35o1w6OHnVPrEXwQIMHGHXnKKZtn8bH2z+2+7iYOTEMjx/uvMCEKMbsThha6wFa6/uAj4C2QD4z9bhPiZsPw9HK3gHdt0GlLrDtSdjyCD46W0V3t27www9w6pQZ6Xb/fvfEWgRvdHqDHrV68PT3T7Ph9w3uDkeIYq8gJYzOSql3gVFAHFDNaVEVkXTccwD/stDhW6j/KhyaQ+Mzz0DSkaz7tGljuj5fv25KGrt2uSPSQrP4WJh37zxqlK3BvV/ey9ELnldSEsKVClKHMQsIAdYBW7XW8ue7t1M+cMdY6PAtQSl/mHqNxB+y7tOokRm0MDDQDCeycaMbAi28MoFlWDJwCddTr3P3/Lu5cuOKu0MSotgqyCOp6sA4IAl4SCn1udOiEsVL1V7sqPAxBFWB2O7wy4Ss9Rp16sCmTWbwwq5dzUx+HqRueF0+v/dz4k/EM3TJUJlDQ4hcFKjSG3gXaAfsBcY6KSZRDF31rQrdfoTqA2D3aNhwL9zINJhftWpmIqa6dU0/ja++cl+whdCzdk/e7PQm8/fMZ+Kmie4OR4hiySsrvYWT+IaYZrdN3oHjS2BFS7iw7+b2ihVh7Vpo2dL0CJ8xw32xFsKoNqMYVH8QL65+ke9+/c7d4QhR7HhlpbdwIqXg1hGmo1/yOTNU+rFvbm4vUwZWrDCPph57DCZNcluoBaWUYkafGTSu3JgHvnmAfWf25X+QECWIVHqLwomIMUOKhN1uHk/Fj4a0VLMtOBgWL4b77oN//QteftljRroN9gtm4f0LCfQNpO/8vpy/dt7dIQlRbBSm0vsyUuktAIIjofM6qPU4/N8EUyF+7YzZ5u8Pn38Ojz4Kr78OzzwDaZ4xF0X1sOp8PeBrDv95mEFfDyI1PREKUcIV5JHUQszjqDLAV8CDTopJeBJLALT4CFrOhFMbTNPbczus2yxmOPTnnoOpU+Gvf4UbN/I+XzHRtnpbpvScwvKDy3lx9YvuDkeIYqEgj6QGA5OA88B9wJfOCEh4qJqPmImZ0LCyDRyaY9YrBW+9ZUoZn31mhhW5ds2dkdptWNNhPNnsSd7a/BbzfrYx9LsQJUxBJlBaCBwGNmutX3FSPLlSSt0N/AUoDczUWq90dQwiH+WbmXqNTYPgx7/B2a3QZDJY/OHFF02F+FNPQc+epo6jVCm7Tz15QjwpKSkw3FnB5/K+3Sez5/Qehi4ZSp3ydVz75kIUMwUpYezXWj8G1C/omyilZimlTiml9mRb310ptV8pdVAp9UJe59BaL7K+/98BmS6tuAqsAB2Xw22jzMCFqzrAleNm25NPwqefmv4anTvD2bPujdUOfhY/Fty3gIiQCO754h6SkpM4mHSQE0k2hn8XwssVJGE0V0q9DdRRStVWSqkCHDsH6J55hVLKAkwFegD1gEFKqXpKqQZKqaXZloqZDn3Zepwornx8ofFEaPsVXNgDy5vAqfVm20MPwTffwO7d0KED/PGHe2O1Q4WQCiwauIhzV8/x08mfuJx6mfHrxrs7LCFcriCtpFoA7wEzgIcpQB2G1no9cC7b6hbAQa31Ia11MjAf6Ku1/llr3SvbckoZE4HvtdY77X1v4UbV+0O3LeBXBlbfBfsmm+a1ffrA99/D77+bkW4PHXJ3pPlqVKkR73Z9lxvW2Qhn7JrBgbMH3ByVEK6l8hs3RykVp7Vunel1KaCW1rpAQ5MqpaKApVrr+tbX/YHuWutHra8fBlpqrZ/O5fhngb9iepjHa60/ymW/YcAwgIiIiKbz588vSJgZkpKSCA0NLdSxjtJoWGcA4qevcmscULTrYUlL4rbzEwi/tomTQZ3YHzaSNJ8gSu3dyx0vvECanx+7336bK9HRuZ4j9NFeaK25PNN9PbDf/fVdliQuybKudmht6ofVp0HpBjQIa0B4QLjL4ikO31GJo/jF4Ig4OnbsuENr3SzHBq11nguw0/rvO5nWxeV3nI3zRAF7Mr3uD8zI9PphYEpBz5vX0rRpU11Ya9euLfSxjnK+tkX/Wdvi7jC01g64HmmpWu95Xeu5SuvvGmh98aBZ//PPWleurHW5clpv2ZLr4ZduteiLdX2KFkMR/HHxDx34WqBmLBmL7zhf3WZmGx38enDGuhrv1dCDFw7W07dP13tP79VpaWlOi6k4fEe1ljiKWwxaFz0OYLu2cU+1p5WUUkpFYDrrjbSeLKjQqeum42QdXiTSuq7IlFK9gd61atVyxOmEIygfuP1FKNsUNj8Ay5vBnXOhfk8zJHrnztCpEyxZYqZ+LWbGrx9Pms7a8dBH+dAwoiFr/7qW+BPxbDi6gY1HN/L9ge/53+7/AVA+qDxtq7elbfW2tKvejsaVG+Nv8XfHRxCiyOxJGKOBDcA84F2l1K84Zi7wbUBtpVQ0JlEMBB5wwHnRWn8LfNusWbPHHHE+4UBVukH37bChH6zrBQ1ehfpjTNLo2hV69IAvvzT1HMVIXEIcyanJWdYlpyazOWEzfhY/mldtTvOqzfln63+itebAuQNsPLoxI4ks3r8YgCDfIFpGtqRttba0u6UdrSJbUTqgtDs+khAFlm/C0FovB+oAKKVaYzrtDS3Im1iHEYkBwpVSCcCrWuuZSqmngRWABZiltf6lYOELjxQaDV02w9bH4eexcHY73PkprFtnEka/fjBnjmlRVUzsetxU2cXMieH8+fPED4/PdV+lFHXK16FO+To80vgRAE4knWDT0U0ZCeSNjW+QtiEto5SSXgppW70tVUpVccVHEqLACtJxD611HGak2gLRWg/KZf0yYFlBz5cfeSTlAXyDoPV/Ibwl7BhuHlG1XwirV0PfvvDww3Dhguno5wUqhVbi3nr3cm+9ewG4dP0SW45vYcPvG9h4bCMzd83kg60fAFCjbA2TPKylkLrl61KwVuxCOEeBEoankEdSHkIpqPMUlG0EG++DFa3MmFTLlsH998PTT5ukMXq0uyN1uFIBpehcozOda5iWcDdSb0g9iCj2vDJhCA9ToY0ZUmTjANg8COqOgC8/h6HD4KWX4M8/UTfSCPxDw4kTUKmSuyN2uMLUg1RLq0ZytWSpBxEu45UJQx5JeaCgytBpDex8Dva/C3/uhI8/N+NPTZpEoD/4JAPjx5uRb72cPfUgcxPn8ulnn0o9iHAZr0wY8kjKQ/n4QbP3oHxz2DoMVjaHV78CiwXL+++bfWbNgjFjvLKUkZ/s9SDLVi3Dv4a/1IMIl/HKhCE8XPRDUKY+rO8HqzvA6dZoQIEZGv2ee2DNGghyRHcgzxXsG0xMjRipBxEu45UJQx5JeYGyjUx/jW/vhQWxZPnb+McfoUYNeO01MymTr1d+jQusMPUg7aq3o231tlIPIuzilf/T5JGUlwgoB6vrAuuBTL2sfX3NdK+PPgrvvANvvgm9e5tWV07mrnk5CsOeepDXN7xOms7aHyQ9iVQuVdnNn0AUN16ZMIQX+XEL3Mg2F3hKClSpAh9+aCZm6tsX2rSBiRPNvyJXtvqD/JjwIxuPbsy1HiQ9gUg9iJCEIYq3zctgSQ1IzTat6y23Q8seZgiRWbNg7FgzVHrfvvDGG1CvnlvC9TSlAkrRpWYXutTsAph6kF0ndpkEkq0eJDw4nDbV2mQkEKkHKXm8MmFIHYYX+Xk8ZBv0Dyzw+1w4Hw9tvoDHHzfDiEyebEoZDRrA3/5mkkhkpOtj9mB+Fj9aVG1Bi6otstSDpLfEkv4gJZsjBhEsdrTW32qth4WFhbk7lCKxWMDik/d8JV7vbBykJWdbmQqhNeH6aVjRHA5+AsHBppPfoUPw7LPwv/9B7drwwgtw/rw7IvcK6fUgQ5sMZXbf2Rx45gCJIxP56r6vGNZ0GJeuX2Lu0bl0+6wbZSeWpcnHTXj2+2f56pevSLyU6O7whYN5ZQlDeJEeZtC/pNt80VpTal/qzW1XT0Dcw6bPxonV0OJjCA+Hd981SeOVV+Ctt2D6dJNMnnoKAgPd9EG8R6XQSvSv15/+9foDpj+IX7Sf1IOUAJIwhOcKqgQdV8D/vQU/vQzntkGb+abjX3Q0fPopjBxpxqJ67jl47z3TU/yhh0zxTThEsG8wMTVjpB6kBJCEITyb8oHbX4CK7WHTIFh5JzSaALeOMNsaNTLzh69ZA6NGwZAhMGkSTJgAPXu6pCluSVOYehDpD+IZvDJhSKV3CVThTugZDz8OhV3Pwck10GoOBFYw2++6C7ZuhQULTFPcXr2gfXtTSd6qlTsj93qZ+4MMbWKm0jmRdCKjBCL9QTyHVHoL7+FfFtp9Dc2mmjqN7xvCybU3t/v4wIABsHevGcBw3z5o3Rr694f9+90XdwmUXg8yuftktg/bzvnnz7PyoZW83O5lygaVZeaumQxYMIAq71Sh5vs1+euivzJj5wz2ndmHmSVauINXljBECaYU1HnSDJm+6X5Y3clMAVt/DPhYv+5+fvDkkzB4sOkp/vbbsGiR6Tn+6qtQWf6idTXpD+IZJGEI71S2oZljY/vTsGecKWm0mQfBmfplhIaallR//7upDP/oI1NRPmIE/OtfICVUt7FVD/Lr2V8zWmJJfxD3kIQhvJdvCLSaDRGdYNsTsKyhqdeI7J11v4oV4YMPYPhwePlleP11kzxefhmeeAICAtwRvchEKUXd8LrUDa+bUQ+SeCmRTcc2ZZRCss8Pkl4CkXoQx/HKOgwhsoh+CLrvhJAoWN/HzCGeej3nfjVrwuefw/bt0LixKWnceit89pkZ7BDwTUmj7vGrZuY/4VaVS1XOUg/ybZtvs9SDzNg1I0s9yJBFQ6QepIi8soQhraREDqVrQ9fNEP887H8PTm0wfTZK1865b9Om8MMPZnn+eXj4YdMUd+JEKp25Rui1tBIz858nya8/yLIDy/jv7v8Cph4kfYIpqQexn1cmDBneXNhkCYCmkyHiLvjxb7C8CTT/CKIftL1/ly7QqRN88YXpKd69O+WxTuQ0a5bp13HLLa6LXxSIPfUgi/YtAkw9SKvIVhmPsDy5HiRmTgznz58nPibe4ef2yoQhRJ4i+0DP3bDpAYh7CE6ugmZTTJ1Hdj4+MGgQ3Huv6cuxaZNZf+0aREWZoUiqVs17KVdOOggWA/bUg2TvDyL1IFmVuIRx48YNEhISuHbtWp77hYWFsXfvXhdFZVvalKUA+Lg5Dsh6PQIDA4mMjMTPz8/NURVBcCR0WgN7xpvlTBy0+RLK3mF7/7NnYceOrDP/+fpCjx7w559w/Lip+zh1KuexgYFm/o68kkqVKuAvj0RcLb0eJH1crOzzg8zYNYP3t5r55GuUrZElgZTEcbFKXMJISEigVKlSREVF5fnLvnTpEqVKlXJhZDmlpl4GwHLbbW6NA25eD601Z8+eJSEhgejoaHeHVTQ+vnDHvyEiBjY/CCtaQNN3odbfc5YIxo/PqPi+ebwPlCplRsZNl5wMiYkmgdhatm0zfT5s/cFSoULOJCKlFZfKrz9IfvUg3q7EJYxr167lmyxE7pRSlC9fntOnT7s7FMeJ6Ag9dkPcX2Hbk6aXeMtPTM/xdHFxJhlklpwMmzdnXefvb+o18qrb0PpmqSS3ZetWsHWNM5VWbvPzg6VLpbTiRAWtB7k19FZ66V4eXw+SmxKXMABJFkXkldcvsALELIV970L8C3BuO9z5OVRobbbvymOY9YJSypQUypUzkz3l5vr1PEsrpfftM3Uq1200Ec5eWrG1lC0rpZUCyq8e5PtfvvfqepASmTCEsEn5wG0joUI72DQQVrWDO16DeqPMNlcLCDAV61FRNjdviY0lpkMHOHcO/vij4KWVoKD861YqV5bSSj4y14PcHXg3TVs39dp6EEkYQmQX3sJM3LR1GOwebUa+bf0pBEW4O7KclILy5c1ShNIKW7aYf22VVipWzDOp+F66ZB6zedCNz5mKWg9SnPuDeGXCKO4d95RS/POf/+Q///kPAJMmTSIpKYmxY8fa3P/06dP06tWL5ORk3n//fdq1a2f3e82dO5eJEyeaxyilSjFt2jQaNmzIkSNH6NWrF3v27HHER/I+/mGmY1+lLrDjWTPybetPUaGawHvSzGx/QZXcHaX98imtAOamf+5c7knl2DH48Uc4cybLYW3BlFbyqqwvwaWV/OpBNvy+wWP6g3hlwijuHfcCAgL45ptvGD16NOHh4fnuv3r1aho0aMCMGTMK/F7R0dGsW7eOsmXL8v333zNs2DC2bNlSmLBLHqWg1qMQ3tqMfLu2G4EDNT4RmKa4zb2sp3fm0soduTQvBlMKyfQI7OC6ddQKCrqZWH780WzPXlpRyr66lTJlvLq0kl89yIajG4ptPYhXJgx7DR8O8fG2t6WmBhVqFs9GjWDy5Lz38fX1ZdiwYbz77ru8/vrree4bv28/o0aN5urVq2zfvp24uDiCgoLsjufOO+/M+LlVq1YkJCTk2OfQoUPce++9TJ8+nebNm9t97hKjzO3QbStsfRwLn5l1Bz8GNITdDiHREBoNIbeAb7BbQ3WJgAAzBa61WXVCRAS1YmKy7qO16btSwNIKkLW0klfdiif3A8omr/4gG45u4JOdn9hdD3I95ToHkw5yIukElUIdWwou0QnDnZ566inuuOMORo0aled+jW6ty7hx49i+fTtTpkzJsX3EiBGsXbs2x/qBAwfywgsvZFk3c+ZMevTokWXd/v37GThwIHPmzKFhw4aF+CQlhG8w+JZCp4KyADoVDkwHsrWWCoywJpCoTIkkyvwbXB2K8fNph1LK9IIPD4e8vlfXruVdtxIXZ/7N3qRZqRx1K7ckJ8Phw1kTS1iYR5ZWilIPUjV2B1tq3WD8uvFM/YtjS8ElOmHkVRK4dOmqUzvulS5dmsGDB/P+++8XqMSQ3bvvvmvXfmvXrmXmzJls3LgxY93p06fp27cv33zzDfXq1St0DCXC1UQ4PNski3QWf+gaBylX4PJhSDoMl4+Yf89uhaMLQKdkOkBBcNWciST936BI8ClEsdaTBQZmKa3YlF9p5fffYfNmos+eNWN8ZRYcnH9ppVKlYl9aKUg9CLUABbPiZzGmwxiHljJKdMJwt+HDh9OkSRP+9re/Ffoc9pQwfvrpJx599FG+//57ypcvn7FPWFgY1atXZ+PGjZIw8vPzeNDZenrrVDg43dRlpPfXyCwtFa4ev5lE0hPK5cNmQqcrCUCmYbaVLwRXM8kjPZFkTi5BldzTvNfd7CytrF+5kva1auWeWDZtMnUrtkorERH5NzEuRqWV3OpBhiwawqqDK0lTkKbTHF7KkIThRuXKlWPAgAHMnDmTRx55pFDnyK+EcfToUfr168enn35KnTp1smzz9/dn4cKFdOvWjdDQUB544IFCxVAinI2DtGw3mrRkOLPZ9v5gSgsh1c1SsX3O7anJcOWYtXRyJGsp5fh3cC3bnBs+AaaeJDQaQqKpdikNfj91M6EEhBebG5o7pPn7Q40aZsmN1qbeJLekcuSISSznzuU81o7SikotQodOB1h/dD1p1r8pklOTmR0/26GlDEkYbjZy5EibdROOMm7cOM6ePcuTTz4JmAr37du3Z2wPCQlh6dKldOnShdDQUPr06eO0WDxaDwf29E5n8YdSNc1iS8pVuPy7SSTZSynntlPz+lnYNP3m/r4hpkSS/VFX+s/+ZYoes6dLb6lVoYJpoZKbq1dtd4ZMX5dLaaV9emklv8dgpUs7PLmPXz+etGyl4FSd6tBShiQMN0hKSsr4OSIigitXruS5/5AhQxgyZEih3mvGjBk2m+NGRUVl9MEoU6YM27ZtK9T5XeWgDiUlJYVm7g7ElXyDIOxWs9iwYc0y2jWplvNxV9JhOLUOUi5lPcCvTLa6k8yV81G2h3cvqYKCzAyMNXNJ5mAGo8xWt/L7pk1E+fmZ14cPw8aNtksrISH21a342n+LjkuIIzk1mUqXYP4CuL8/nCyVzOaEPErBBSQJQwgPleoTDGUamCU7rSH5z5slk8yPvS7ug8TlkHo16zEBFXImlIyfbzETUImbfHxylFaO1KlDVPYmxrmVVtKXDRvM9hs3cp7f3tIKsOtxUwo+UyaA8heSOXH1SZgkraS81uuvv85XX311c8W1K9zbtTNj3i9Rf1cLR1AKAsqZpVyTnNu1hmunbNefnNsJCQshLdsNLKiKjUQSRWDKSUhLMcPFi5zsLa3kVbfy22+wfr0Z5Ti70NCbFfZly1L+QrKZt2X2bBgzxpRUHER+w8XISy+9xEsvvZTxOnXP9jz2FqIIlDJjYwVFQHirnNvTUuFaYtamwunJ5fQG+H1eRquxVgBfPGwmpcrRZNj6uCuoSsls4WUvHx/Tr6RiRWicx7waV67kXVpZufLmvqmpDp973mMShlLqNuAfQDiwWms9zc0hCeG9fCwmAQRHAjbGLku7YZoFJx1m387l3Fo14GblfOIKuPpHtvP5m46LuTUZDqxYolt42S04GGrVMkt2iYlQo8bNWSGTkx1eynBJwlBKzQJ6Aae01vUzre8OvAdYgBla6wm5nUNrvRf4u1LKB/gfIAlDCHfx8cu4+Z8I9uHWhjFZt6deg8tHb5ZMMrfyOrYQrmcbbt0SbK2Aj7LdU95f5u7Il61ZIR1cynBVCWMOMAVzowdAKWUBpgJdgARgm1JqCSZ5vJnt+Ee01qeUUn2AJ4BPXRG0EKKQLIFQuo5ZbLmRZG3VdSTnY6/Tm+HG+az7+5XOfciVkCjwc+90ysWCvbNCFoHSWue/lyPeSKkoYGl6CUMp1RoYq7XuZn09GkBrnT1Z2DrXd1rrv+SybRgwDCAiIqLp/Pnzs2wPCwvDnmHPU1NTsRRm9EEHCv59PwBXbqnr1jgg5/U4ePAgFy5ccNn7hz7aC601l2d+57L3tKXp3zuhgZ0frXZrHGCaZ4eGhro7DKfE4ZuWRGBqIoEpJ8y/qScITDlBkPVni846J/oNn9JcpiLJ/lW5ZqnENd/K5l9LZa75RpCmXNPCqzj8Thzxf6Vjx447tNY5Wtu4sw6jKnAs0+sEoGVuOyulYoB+QACwLLf9tNbTgekAzZo10zHZmrjt3bvXrjGiLl265LSxpOydDyO9a9i1a9cKPR9Gum3bttG6dWvmz59P//79iY2NZdKkSSxdutSu47Nfj8DAQBrnVTnnYPG+vqSkpJD99+lqSUqB1m6PAyA2NrZkxqE1XD+Tpe+JX9Jh0o7tpKJfIly20Ss/sJKNvifpdSnVzSM2BygOvxNn/l/xmEpvrXUsEOvmMBzClfNhgCkdPP/883Tt2rVQx4ubSmQHwuJGKTMHe2AFMzui1U9XrDdrbZ3gKvuAkJcPw5k4OPqFGQcs43w+EFQ1Z5PhjDG8qpa8QSFz4c6EcRyolul1pHVdkdk7497w5cOJPxFvc1thH0k1qtSIyd0n57mPK+fDAPjggw+49957c+3NvW3bNoYNG8aCBQuomVdbcSE8gfKB4CpmqdAm5/a0FDMoZJKNhHJyNVw5TpZBIX38TAuv7EOupJdWAiuVmAp5dyaMbUBtpVQ0JlEMBBwy+l1xn3EPXDcfxvHjx1m4cCFr1661mTA2b97MM888w+LFi6levXrhP5AQnsLH1/RcD7kFImJybk+9blp4ZR5qJT2pHP8Wrp3Mur8l0PpoK4raF/zh/7bdfOwVEgUB5V2aUGqpJLSvc+qmXdWs9nMgBghXSiUAr2qtZyqlngZWYFpGzdJa/+Kg97OrhJFXScCZdRjguvkwhg8fzsSJE/Hxydlpau/evQwbNoyVK1dSpUqVQscghFexBEDp2maxJeWKNYEcyZFQKl49APFLsu7vG2q770n6Yy+/4jNnd35ckjC01oNyWb+MPCqwi/B+xb6EAa6ZD2P79u0MHDgQgDNnzrBs2TJ8fX0pU6YMlStX5tq1a+zatUsShhD28g2GsHpmyWZTbCwxdza20Vz4yM15UFKSsh7kXzbvJsMFnPZXhWoC77HW4wTJFK1ewxXzYRw+fDjj5yFDhtCrVy/uvvtuYmNjKVOmDDNnzqRLly6EhIS4vXWHEF7BPwz8G0JZG5M9aQ3J53IOCJl0GC78H/yxzHR6zCywou1EEhJtWnhlGxTSr10alurAnvFmci8H8sqEYe8jqeLA2fNh5CciIoKlS5fSo0cPZs2aRcuWubZsFkIUlVKmTiOgPJS30dZOp5k6kvREkrmH/Nlttqf9zTwoZEB5/Bpbh+06NBvqj3FoKcMrE0ZxfyTlyvkwMpszZ07GzzExMRkliurVq/PLLw6pPhJCFIXygaDKZinotL+nYuHKsZvzzutUh5cyvDJhCCGEV8pr2t+ribCkxs1HWmnJDi9lSMIoRmQ+DCFEof08PmPI+QwOLmV4ZcLwpDqMzGQ+DCFEoZ21MSRKWjKckSla81Tc6zCEEMLhepgpWuM/K2OGrxmSlM8BBSdTYAkhhLCLJAwhhBB28cqEoZTqrZSa7sr5GoQQwtt5ZcLQWn+rtR4WFhbm7lBsUkoxcuTIjNeTJk3KMRdGZqdPn6Zly5Y0btyYDRs2FOi9Lly4QO/evWnYsCG33347s2fPBsy4/b169SpU/EKIkskrE0Zxlz4fxpkzZ+zaP30+jF27dhV48qSpU6dSr149du/eTWxsLCNHjiQ5+zSOQgivMTylEY+ed04LUa9sJWW3HcPhz3ibm4JSU6EwU7SWbQRNJ+e5iyvnw1BKcenSJbTWJCUlUa5cOXx9s/7aZT4MIYQ9vDJheEI/DFfNh/H000/Tp08fqlSpwqVLl/jiiy+yDHUu82EIIezllQnD7n4YeZQErnrJfBgrVqygUaNGrFmzht9++40uXbpkPNaS+TCEEAUhdRhuNHz4cGbOnMnly5cLfY4RI0bQqFGjHMuECRMAmD17Nv369UMpRa1atYiOjmbfvn0AVK5cmcDAQHbt2uWQzyOE8G5eWcLwFK6YD6N69eqsXr2adu3acfLkSfbv30+NGjXYs2ePzIchhCgQKWG42ciRI+1uLVUYY8aMYfPmzTRo0IBOnToxceJEwsPDM7anz4fx1FNPsWXLFqfFIYTwfFLCcANXzodRpUoVVq5cmWO9zIchhCgoSRjCIwx/oRHnz58n3t2BCFGCeWXC8IRmtbbIfBhCiOLMKxOGpw5vLvNhCCGKKnZILLGxsU45t1R6CyGEsIskDCGEEHaRhCGEEMIukjCEEELYRRKGG7hyPox9+/bRunVrAgICmDRpUpZty5cvp27dutSqVStjKBGAqKgop3YmFEJ4Jq9sJWW34cMhPt7mpkIPb96oEUyenOcu6fNhjB49Okuv69ykz4cxY8aMAodTrlw53n//fRYtWpRlfWpqKk899RQ//PADkZGRNG/enD59+lCvXr0Cv4cQomSQEoYbZJ4PIz9mPoxRLF68mEaNGnH16tUCvVfFihVp3rw5fn5+WdZv3bqVWrVqUaNGDfz9/Rk4cCCLFy/Oss/Vq1fp0aMHn3zySYHeUwjhnbyyhGF3x708SgLOHt7cVfNh5Ob48eNUq1Yt43VkZGSWsaSSkpIYOHAggwcPZvDgwVy6dMmejyWE8GJemTA8oeOeq+bDKKy+ffsyatQoHnzwQaecXwjheeSRlBu5Yj6M3FStWpVjx45lvE5ISKBq1aoZr9u0acPy5cvRWhc6NiGEd/HKEoancMV8GLlp3rw5Bw4c4PDhw1StWpX58+czb968jO3jxo1j3LhxPPXUU3z44YeFeg8hhHeREoabOXs+jBMnThAZGck777zDa6+9RmRkJBcvXsTX15cpU6bQrVs3brvtNgYMGMDtt9+e5dj33nuPq1ev5lvPIoQoGaSE4QaunA+jUqVKJCQk2NzWs2dPevbsmWP9kSNHMn6ePXs2gFR6CyGkhCGEEMI+UsIoRmQ+DCFEcVYiE4bWGqWUu8PIwVPmw5CWU0KUTCXukVRgYCBnz56Vm14haa05e/YsgYGB7g5FCOFiJa6EERkZSUJCAqdPn85zv2vXrrn9pph20rSe8rHsdWsckPV6BAYGEhkZ6eaIhBCuVuIShp+fH9HR0fnuFxsbS+PGjV0QUe6S+jVAa02pfalujQOKx/UQQriXRz2SUkqFKKW2K6V6uTsWIYQoaVySMJRSs5RSp5RSe7Kt766U2q+UOqiUyn2kvJueB750TpRCCCHy4qpHUnOAKcD/0lcopSzAVKALkABsU0otASzAm9mOfwRoCPwfILWtQgjhBi5JGFrr9UqpqGyrWwAHtdaHAJRS84G+Wus3gRyPnJRSMUAIUA+4qpRaprVOs7HfMGCY9WWSUmp/IcMOB4rDtHPhKFU84nD/9QhXI4rJtZDfSWYSR/GKAYoexy22Vrqz0rsqcCzT6wSgZW47a61fAlBKDQHO2EoW1v2mA9OLGpxSarvW2u095iSO4hWDxCFxFPcYnBmHx7WS0lrPcXcMQghRErmzldRxoFqm15HWdUIIIYohdyaMbUBtpVS0UsofGAgscWM82RX5sZaDSBw3FYcYQOLITuK4qTjEAE6KQ7liiAyl1OdADKYi5iTwqtZ6plKqJzAZ0zJqltb6dacHI4QQolBckjCEEEJ4Po/q6S2EEMJ9SnzCUEpVU0qtVUr9n1LqF6XUP6zryymlflBKHbD+W9ZF8ViUUruUUkutr6OVUlusveG/sNb3ODuGMkqpBUqpfUqpvUqp1u64HkqpEdbfyR6l1OdKqUBXXA9bIxPk9vmV8b41np+UUk2cHMfb1t/LT0qphUqpMpm2jbbGsV8p1c2ZcWTaNlIppZVS4dbXTrkeeYwW8Yz1evyilHor03qXXQulVCOl1I9KqXhlhi5qYV3vzO9Gge5bDotFa12iF6Ay0MT6cyngV0znwLeAF6zrXwAmuiiefwLzgKXW118CA60/fwQ84YIY/gs8av3ZHyjj6uuB6adzGAjKdB2GuOJ6AO2BJsCeTOtsfn6gJ/A9oIBWwBYnx9EV8LX+PDFTHPWA3UAAEA38BlicFYd1fTVgBfA7EO7M65HLtegIrAICrK8ruuNaACuBHpk+f6wLvhsFum85KhaH/kfzhgVYjBmuZD9QOdMvZ78L3jsSWA3cBSy1/nLPZLpBtAZWODmGMMyNWmVb79Lrwc2OneUw/YWWAt1cdT2AqGw3BZufH/gYGGRrP2fEkW3bPcBc68+jgdGZtq0AWjszDmABZsieI9xMGE67HjZ+J18CnW3s59JrYT3//dafBwHzXPHdyBZTnvctR8VS4h9JZabM8CWNgS1AhNY60brpBBDhghAmA6OA9F7s5YHzWusU6+sEzI3UmaKB08Bs66OxGUqpEFx8PbTWx4FJwFEgEbgA7MD11yNdbp/f1ogFrorpEcxfjS6PQynVFziutd6dbZMr46gDtLM+olynlGruhhgAhgNvK6WOYb6zo10Zh533LYfEIgnDSikVCnwNDNdaX8y8TZuU7NTmZMoM2X5Ka73Dme9jB19MkXua1roxcBlTtM3goutRFuiLSWBVMOOIdXfme9rLFZ8/P0qpl4AUYK4b3jsYeBF4xdXvnY0vpgTaCvgX8KVSbpl7+QlghNa6GjACmOmqN3b1fUsSBqCU8sNc9Lla62+sq08qpSpbt1cGTjk5jDZAH6XUEWA+5rHUe0AZpVT6EC6u6A2fACRorbdYXy/AJBBXX4/OwGGt9Wmt9Q3gG8w1cvX1SJfb53f5iAXKjKfWC3jQelNwdRw1MYl8t/X7GgnsVEpVcnEcCcA32tiKKZmHuzgGgL9ivp8AX2EGVsXZcRTwvuWQWEp8wrD+RTIT2Ku1fifTpiWYLwLWfxc7Mw6t9WitdaTWOgrT632N1vpBYC3Q34VxnACOKaXqWld1wgwr79LrgXkU1UopFWz9HaXH4dLrkUlun38JMNjaCqUVcCHTIwGHU0p1xzy27KO1vpItvoFKqQClVDRQG9jqjBi01j9rrStqraOs39cETAXsCVx7PRZhKr5RStXBNNA4gwuvhdUfQAfrz3cBB6w/O+1aFOK+5ZhYnFEB40kL0BZTbPsJiLcuPTH1B6sxv/xVQDkXxhTDzVZSNTBf9oOYv14CXPD+jYDt1muyCCjrjusB/BvYB+wBPsW0enH69QA+x9Sb3MDcDIfm9vkxDROmYlri/Aw0c3IcBzHPotO/qx9l2v8laxz7sbbacVYc2bYf4Walt1OuRy7Xwh/4zPr92Anc5Y5rYb2H7MC0zNoCNHXBd6NA9y1HxSI9vYUQQtilxD+SEkIIYR9JGEIIIewiCUMIIYRdJGEIIYSwiyQMIYQQdpGEIYQQwi6SMIRwAuuQ3//J9Po5pdRYN4YkRJFJwhDCOa4D/dLniRDCG0jCEMI5UoDpmMHo8qWUmqOU6p/pdZKzAhOisCRhCOE8U4EHlVJh7g5ECEeQhCGEk2gz3PT/gGfdHYsQjiAJQwjnmowZoC4kn/1SsP5/VEr5YAbWE6JYkYQhhBNprc9hphIdms+uR4Cm1p/7AH5ODEuIQpGEIYTz/QczsU9ePgE6KKV2Y+Yqv+z0qIQoIBneXAghhF2khCGEEMIuvvnvIoRwFKXUS8B92VZ/pbV+3R3xCFEQ8khKCCGEXeSRlBBCCLtIwhBCCGEXSRhCCCHsIglDCCGEXf4fxj9EhEcRSFgAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots()\n",
    "\n",
    "color = ['blue', 'green', 'orange', 'red']\n",
    "\n",
    "for i in range(len(N_f)):\n",
    "    plt.errorbar(N_u, mean_2d[i], yerr=std_2d[i], linestyle='None', marker='^',\n",
    "                 color = color [i])\n",
    "        \n",
    "    plt.plot(N_u,mean_2d[i], label = 'N_f =' + str(int(N_f[i]/1000)) + 'k', color = color [i])\n",
    "\n",
    "ax.set_xticks(np.linspace(20,200,num=10))   \n",
    "ax.set_yscale('log')\n",
    "ax.set_ylim([1e-4,2])\n",
    "plt.xlabel('N_u')\n",
    "plt.ylabel(r'$\\varepsilon_{PINN}$')\n",
    "plt.legend()\n",
    "plt.grid()\n",
    "# fig.savefig('Trend.eps', format = 'eps')"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
